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ABSTRACT 

We perform a cosmological simulation in order to model the growth and evolution 
of Population III (Pop III) stellar systems in a range of host minihalo environments. 
A Pop III multiple system forms in each of the ten minihaloes, and the overall mass 
function is top-heavy compared to the currently observed initial mass function in the 
Milky Way. Using a sink particle to represent each growing protostar, we examine the 
binary characteristics of the multiple systems, resolving orbits on scales as small as 20 
AU. Wc find a binary fraction of ~ 36%, with semi-major axes as large as 3000 AU. 
The distribution of orbital periods is slightly peaked at < 900 yr, while the distribution 
of mass ratios is relatively flat. Of all sink particles formed within the ten minihaloes, 
~ 50% are lost to mergers with larger sinks, and ~ 50% of the remaining sinks are 
ejected from their star-forming disks. The large binary fraction may have important 
implications for Pop III evolution and nucleosynthesis, as well as the final fate of the 
first stars. 

Key words: stars: formation - Population III - galaxies: formation - cosmology: 
theory - first stars - early Universe 



1 INTRODUCTION 

The period during which the first stars arose marked a piv- 
otal turning point in the evolution of the Universe. The first 
stars, also known as Population III (Pop III), are believed 
to have formed at z > 20 within dark matter (DM) mini- 
haloes of 10 6 M 



z 



. .Haiman et al.lll996l : | 
19971 : lYoshida et alj|2003h . Afterwards they generated ion- 
izing photons that contributed to the reionization of 
the Universe (e.g. Kitavama et alj 2004 : [Sokasian et al 



2007) 



2004IWhalen et al.ll2004 l Alvarez et al .1120061 :1 Johnson et al 



while those Pop 
lives as supernovae (SNe) 



III stars that 
contributed to 



ended their 
the metallic- 



Mori et al. 2002 


; Bromm et al.1 2003: Wada & Venkatesan 


2003|; 


Norman et all 2004 


:lTornatore et al.ll2007l: iGreif et al. 


20071, 


20ld;IWise & Abel 2008: Maio et al. 201 1|; recently re- 


viewed inlKarlsson et al. 


2012). This metal-enrichment and 



growing radiation background provided by the first stars de- 
termined the environment in which later stellar generations 
formed. 

The mass of a given Pop III star plays the central role in 
determining how the star will influence its surroundings. For 
instance, more massive stars are more efficient in produc- 
ing ionizing photons, per unit stellar mass, than low-mass 



E-mail: athcna.stacy@nasa.gov 



stars l|Bromm et al j200ll : ISchaererl2002l h Furthermore, non- 
rotating primordial stars with mass 140 Mq < M„ < 260 Mg 
are belie ved to have exploded as pair-instability supernovae 
(PISNe: iHeeer fc Wooslevl 12002? ), releasing the entirety of 
their metal content into the IGM and surrounding haloes, 
while stars within the range 15 Mq < M» < 40 Mq ended 
their lives as core-collapse SNe. On the other hand, non- 
rotating Pop III stars with main sequence masses in the 
range 40 M Q < M, < 140 M Q or M, > 260 M Q are expected 
to collapse directly into black holes (BHs), thus contributing 
no metals to their surroundings. 



BH remnants of Pop III stars may emit copious amounts 
of radiation during subsequent mass accretion, however. 
While this radiation will heat and ionize its local surround- 
ings, the X-ray component has a long mean-free path that 
allows it to ionize the IGM and minihalo gas as far as 
10 — 100 kpc away. This contributed to reionization, while 
the boosted electron fraction can catalyze the formation 
of additional H2 molecules, mildly e nhancing the globa l 
cooling and star-format i on rate (e.g., Madau et al. 2004 ■ 
i Ricotti fc Ostrikerl 12004 iKuhlen fc Madaul 120051 : iHaimanl 
l201ll : iJeon et al.ll2012? ) . Recent work has demonstrated that 
a BH remnant which has a stellar binary companion, giv- 
ing rise to a high-mass X-ray binary (HMXB), will exert 
stron ger feedback on its surroundings tha n an isolated BH 
(e.g- lkirabel et al.|[2uTll : |jeon et al.ll2012t ). This is because 
accretion onto solitary BHs is more susceptible to feed- 
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back that in turn reduces the BHs accretion rate and lu- 
minosity (e.g. Johnson fc Brominll2007l ; lAIvarez et al.ll2009l ; 



iMilosavlievic et al.ll2009l ) , whereas the near-Eddington lumi- 
nosity of a HMXB can persist through Roche lobe overflow. 
In this way the binary characteristics of Pop III stars will 
influence the nature of Pop III feedback. Also intriguing are 
observations indicating that ultra-luminous X-ray sources 
(ULXs), which may be accreting, intermediate-mass BHs, 
occur at a hi gher rate in less m assive and lower-metallicity 
galaxies (e.g. lSwartz et al1l2008f ). Recent modeling of binary 
system evolution over a range of metallicity also support a 
high er rate of ULX form ation in lower-metallicity star clus- 
ters (I Linden et al.|[201ol ). We note that BH binaries are de- 
tectable in c hannels other than el ectromagnetic radiation. 
For instance. iKowalska et al.l (|2012T ) find that if Pop III stars 
left behind a significant BH population, then their gravita- 
tional wave (GW) emission will dominate the low-frequency 
(< 100Hz) spectrum as long as the Pop III binary fraction 
was at least 10" 2 . 

A binary companion can also influence a star's evolu- 
tion during and after its main sequence (MS) stage, as well as 



the t y pe of death the star w ill under go (e.g. Wellstein et al.l 
200ll: IPetrovic et al l [iooBI 2005b: ICantiello et alj 120071 : 



de Mink et al.ll2009l ; lLanger!l2012i ; ISana et alj|2012h . For in 



stance, massive stars that leave behind BH remn ants may 
power collapsar gamma-ray bu rsts (GRBs; e.g. IWooslevI 
ll993l ; lLee fc Ramirez-RuizfeOOd ). However, this requires suf- 
ficient angular momentum in the Pop III progenitor for 
an accretion torus to form around the remnant BH. The 
progenitor star must also lose its hydrogen envelope to 
enable the r elativistic jet to p enetrate through and exit 
the star (e.g. IZhang et al.ll2004 but see lSuwa fc Iokall201ll ; 
iNakauchi et alj|2012f ). One way to fulfill both of these condi- 
tions is for the GRB progenitor to have a close binary com- 
panion that removes the hydrogen envelope due to the heat- 
ing during a com mon-envelope phase (e.g. iLee et al.l [20021 ; 
llzzard et aj]l2004l ). 

A star which receives mass from or merges with its bi- 
nary companion potentially undergoes a process of 'rejuve- 
natio n' in which its central hydrogen abundance increases 
(e.g. iHellingsl Il983l ; H rami fc Langerl 1 19951 ; iDrav fc Tout] 
l2007h . This causes the star's apparent age to be less than 
its true age, and may lead the star to app ear as a blu e 
straggler on the Hertzsprung-Russel diagram (|Langerll2012l) . 
Whether rejuvenation will result from mass transfer de- 
pends on the evolutionary stage of the binary components, 
and is not expected to occur if one of the members is al- 
ready post MS. For instance, a star wit h an undermas- 
sive helium core may result instead (e.g. iBraun fc Langerl 
Il995l ). Stellar mergers may furthermore lead to the gen- 
eration of large-s cale magnetic fields observed within mas- 
sive M S stars (e.g. lDonati fc Landstreetll2009l ; lFerrario et al.l 
l2009h . In addition, while single stars must typically have 
masses > 20 M q to g o through a Wolf-Rayet (WR) phase 



(|Hamann et alj 20061). and even lar ger masses at lower 
metallicities ( Mevnet fc Maederll2005h . a member of a close 



binar y may become a WR star at a lower ma ss of ~ 10 Mq 
(e.g.. [Vanbeveren et al.ll 199 i lCrowtherll2007l ). 

A star that gains mass from its binary companion 
through an accretion disk may spin up to critical rota- 
tion speed, and rapid rotation ca n alter a star's ev o lution 
in a number of ways (see, e.g., IStacv et alj 1201 ll . |2012| ; 



iMaeder fc Mevnet|[2012l ; lYoon et alj|2012h . Rotation will al- 
low for rotationally induced mixing and possibly chemically 
homogeneous evolution (CHE). CHE, in turn, may lower the 
minimum mass at which a star will undergo a PISN to as low 
as ~ 65 M Q (e.g.. lChatzopoulos fc Wheeled 12012] ). If a ro- 
tating star instead leaves behind a BH remnant, a collapsar 
GRB may result even without a binary companion. This is 
because CHE may allow stars to bypass the red giant phase 
to become a WR star. This evolutionary path may further- 
more allow the star to retain enough angular momentum to 
become a GRB (e.g. lYoon fc La ngcr 2005; IWooslev fc Hegerl 
l2006h . Rotationally induced chemical mixing will generally 
enhance a star's luminosity, temperature, and metal pro- 
duction, though this depends sensitively o n the assumed 
magnetic field evolution within the star (e.g., lEkstrom et al.l 
120081 ; lYoon et alj l2012i ) . Another possible consequence of 
both rotation and mass transfer is enhancement in surface 
nitrogen abundance (e.g. lEkstrorn ct aT, 200§; IL"anger et all 
120081 ; iBrott et alj l201ll ; lLangerll2012i ). 

Observations of massive stars in the Milky Way find 
that ~ 45-75% of O-stars have spectroscopic b inary com- 
panions ( Mason et al 1 l2009l ; ISana fc Evansll201lJ) , and calcu- 
lations of binary evolution furthermore fi nd that 20-30% o f 
O-stars will merge with their companions (|Sana et al.ll2012T l. 
Binarity therefore plays a crucial role in the evolution of 
massive stars in the Milky Way. While similar observations 
cannot be made for Pop III stars, we may still study them 
numerically. Early computational studies of Pop III star for- 
mation indicated tha t they would grow to be highly mas- 
sive (> 100 Mq : e.g. I Abel et a.1.1 120021: iBromm et alj |2002| ; 



sive iuu iviq| e.g. lApei ex auiz uuz; rjromm ex amzuuzi ; 
IBromm fc Loebl 12004 lYoshida et all 120081 ). and that each 



minihalo would host only a single Pop III star. More recent 
work, however, has found that primordial star-forming gas 
will undergo continued fragmentation well after the initial 
star has formed. Pop III multiplicity has been found both 
in small-scal e simulations that u tilize idealized initial con- 
ditions (e.g. iMachida et aljlioosl , Clark et al 2008; 2011a), 
as well as s imulations that are initialized on cosmological 
scales, (e.g. iTurk et al.l|200d; IStacv et alj|20ld : Clark et al. 



2011b. iGreif et al.ll2011 



2012!) and also when accounting fo r 



radiative feedback (e.g. lSmith et alj|201ll ; IStacv et ai1l2012l ) . 
Recent numerical advances have thus shown that the pro- 
cesses that lead to bin ary formation apply to both enriched 
and metal-free gas (e.g. lTohlinell2002j ). and a new picture has 
emerged in which Pop III stars typically form in multiples 
and within disks. 

In order to determine the role of multiplicity in Pop III 
evolution and feedback, in this study we aim to better con- 
strain the nature of Pop III binaries. We do this through 
simulating a large cosmological box which contains approx- 
imately ten minihaloes in which primordial stars form by 
z ~ 20. We examine the rate of binary formation and pro- 
tostellar mergers, the typical characteristics of the binaries, 
and the overall mass function at this stage of evolution of 
the Pop III systems. Stars ejected from their Pop III cluster 
may stop accreting and remain low-mass and long-lived to 
the present day. Such stars may be observed in dwarf galax- 
ies or the Milky Way halo, so we furthermore discuss the fre- 
quency of low-mass sink ejections in our suite of minihaloes. 
In Section 2 we describe our numerical methods, while we 
present our results in Section 3. We address our caveats in 
Section 4, and conclude in Section 5. 



© 2012 RAS, MNRAS 000.1111141 



Pop III Binaries 3 




Figure 1. Temperature versus number density of the gas just prior to the formation of the first sink particle in each minihalo. The 
similarity in the early thermal evolution is striking. At n > 10 s cm -3 , three-body reactions rapidly convert the gas into fully molecular 
form. The corresponding boost in H2 cooling balances adiabatic heating, resulting in nearly isothermal evolution. H2 line opacity will 
become important only at even higher densities, eventually causing deviations from the near-isothermality seen here. 



2 NUMERICAL METHODOLOGY 
2.1 Initial Setup 

We perform our investigation using GADGET-2, a widely- 
teste d three-dimensional N-body and SPH code (jSpringell 
120051 ). We begin with a 1.4 Mpc (comoving) box contain- 
ing 512 3 SPH gas particles and the same number of DM 
particles. The simulation is initialized at z = 100. Posi- 
tions and velocities are assigned to the particles in accor- 
dance with a ACDM cosmology with Qa — 0.7, Qm — 0.3, 
Qb = 0.04, as = 0.9, and h — 0.7. Each gas particle has a 
mass of m sp h = 120 Mq, while DM particles have a mass 
of mDM = 770 Mq. Once a gas particle reaches a threshold 
density of 10 3 cm -3 , we employ 'marker sinks' instead of 
following these regions to increasingly higher densities. This 
allows us to efficiently evolve the entire box through a time 
period over which ten minihaloes form. Though the marker 
sinks are numerically similar to the protostellar sink parti- 
cles to be described in Section 2.4, note that the marker sinks 
are implemented at much lower densities than the protostel- 
lar sink particles of the final highest-resolution calculations. 



The location of the marker sinks coincides with the cen- 
ters of the minihaloes. In particular, we locate the first ten 
sink particles which form within the first ten minihaloes. 
Once these positions are ascertained, the cosmological box 
is reinitialized at z = 100 for each individual minihalo, but 
with 64 'child' particles added around a 100-140 kpc region 
where the target halo will form. Larger regions of refinement 
are used for minihaloes whose mass originates from a larger 
area of the cosmological box. Particles at progressively larger 
distances from the minihalo are given increasingly larger 
masses, such that in the refined initial conditions there is 
a total of ~ 10 7 particles. The most resolved particles are of 
mass m sp h =1.85 Mq and moM = 12 Mq. 



2.2 Cut-Out and Refinement Technique 

The refined simulations are run until the gas reaches a den- 
sity of 5 x 10 7 cm -3 . All particles beyond 10 physical pc 
are then removed from the simulation. Once the gas has 
reached these densities, this central star-forming region is 
gravitationally bound, and the effects of the outer minihalo 
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Figure 2. Velocity structure of the gas just prior to the formation of the first sink particle in each minihalo. Minihaloes are shown in 
order going from top left to bottom right. Solid lines show radial velocity v [lL( j L , dashed lines depict rotational velocity v TO t, and dash-dot 
lines show the free-fall velocity vg. Blue dotted lines represent M tur t, and follow the scaling shown on the right-hand y-axis. Turbulence 
is generally sonic (M tur ]- ) =1), though in some cases the level of turbulence can become mildly supersonic. Over these distance scales, 
u ra< j and v TO t both remain on the order of ~ 1/2 of Vff. While the H2 chemistry leads to very similar thermal histories for each minihalo, 
differences in mass inflow history yield a greater variety in the velocity structure of the inner gas. In particular, the relative magnitude 
of i; rat j and v ro t differs in each case. Though they remain within a factor of two of each other, t> ra[ j dominates in some cases while Vrot 
dominates in others. 



and neighboring haloes can be safely ignored. Also note that 
the gas at the cut-out edge has a free-fall time of ~ 10 7 yr 
and will undergo little evolution over the next 5000 yr fol- 
lowed in the simulation. Our cut-out technique leads to the 
propagation of a rarefaction wave starting from the cut-out 
edge due to the vacuum boundary condition. However, the 
wave will only travel a distance of c 3 t, where c s is the gas 
sound-speed (~ 2 kms -1 ), and the time t is 5000 yr. This 
corresponds to an insignificant distance of ~ 10~ 2 pc (2000 
AU) from the cut-out edge, over two orders of magnitude 
smaller than the 10 pc box size. 

In addition, each remaining SPH particle is replaced 
with 256 child particles, each of which is placed randomly 
within the smoothing kernel of the parent particle. The mass 
of the parent particle is then evenly divided amongst the 
child particles. Each of these particles inherits the same 
chemical abunda nces, velocity, and en tropy as the parent 
particle (see, e.g.. lBromm fc Loebll2003l ; Clark et al. 2011b). 



This ensures conservation of mass, internal energy, and lin- 
ear momentum. After this process, each SPH particle in the 
new cut-out simulation has a mass ofm sph = 7.2 xl0~ 3 M Q , 
and the new resolution mass is M res ~ 0.4 Ma, the approx - 
imate mass contained in one SPH kernel (jBate et al.lll995l l. 
The cut-out and refinement is performed for the central star- 
forming region within each of the ten identified minihaloes. 
We then evolve the star-forming disks in these ten different 
minihaloes over 5000 yr of protostellar accretion. 



2.3 Chemistry, Heating, and Cooling 

We utilize the same ch emistry and therm al network as 
described in detail by iGreif et all (|2009h and used in 
IStacv et all |2012T l. Briefly, the code follows the abundance 
evolution of H, H+, H", H 2 , H+, He, He+, He ++ , and e", 
as well as the three deuterium species D, D + , and HD. All 
relevant cooling mechanisms, including H2 collisions with H 
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Figure 3. Angular momentum profile of gas within each mini- 
halo, measured just prior to the formation of the first sink. The 
thick lines of various colors and linestylcs represent eight different 
minihalocs from our suite. The thin yellow lines denote angular 
momentum profiles found in separate cosmological simulations. 
Solid yellow line is taken from Stacy et al. (2010), dotted yel- 
low line from Yoshida et al. (2006), and dashed yellow line from 
Abel et al. (2002). Thick black dashed line shows an approximate 
powerlaw fit to these profiles, j oc M enc Note the near conver- 
gence of each of these angular momentum profiles, including pro- 
files obtained from varying cosmological realizations and differing 
techniques for hydrodynamics computation. 



and He as well as other H2 molecules, are included. The 
thermal network also includes cooling through H2 collisions 
with protons and electrons, H and He collisional excitation 
and ionization, recombination, bremsstrahlung, and inverse 
Compton scattering off the cosmic microwave background. 

Further H2 processes are also included to properly 
model gas evolution to high densities. For instance, the 
chemistry and thermal network includes three-body H2 for- 
mation and the concomitant H2 formation heating, which 
become important at n > 10 8 cm -3 . T hree-body H2 for - 
mation rates are uncertain, however, and iTurk et all (|201ll ) 
determined that the variation in suggested rates leads to sig- 
nificant differences in the gas collapse at high density, includ- 
ing long-term disk stability and fragmentation. For our sim- 
ulations, we choose the intermediate rate from Palla et al.l 
j 19831 ). When n > 10 9 cm" 3 , cooling through H 2 ro- 
vibrational lines becomes less effective as these lines grow 
optically thick. We account for this using an escape probabil- 
ity formalism together with the Sobo lev approximation (see 
lYoshida et aHl200fJ : iGreif et al.ll201ll for further details). 



2.4 Sink Particle Method 

SPH particles that reach densities of n max = 10 13 cm -3 are 
converted into sink particles. Sinks grow in mass by accret- 
ing surrounding gas particles that lie within a distance of 



where r. 



0.5 



Mr 



1/3 



20 AU, and 



Along with the criterion that the accreted particles 
must be within a distance d < r acc from the sink, we 
also require that the particles are not rotationally sup- 
ported against infall onto the sink. This corresponds to 
Jsph < Jccnt, where j'sph = v Iot d is the specific angular mo- 
mentum of the gas particle, j CC nt = V GM s i n kr aC c is the an- 
gular momentum required for centrifugal support, and w ro t 
and d are the rotational velocity and distance of the particle 
relative to the sink. In addition, a sink accretes all particles 
which come within d < r m i n = 4AU. These same criteria 
are used not only for gas particles but also neighboring sink 
particles, so our algorithm allows for the merging of two sink 
particles. When the sink first forms, it immediately accretes 
the majority of particles within r- acc , giving it an initial mass 
of M rcs — 0.4 Mq. Each time the sink grows, its position and 
velocity are set to the mass- weighted of that of the sink and 
the accreted particles. 

After a particle becomes a sink, its density, temper- 
ature, and chemical abundances are no longer updated, 
though its position and velocity continue to evolve through 
gravitational interactions. Each sink is instead held at a con- 
stant density n max and temperature of 1000 K. The sink thus 
exerts a pressure on the surrounding particles, avoiding the 
formation of an artificial pressure vacuum around its a ccre- 
tion radius (see lBromm et al.|[20o3 ; iMartel et al.ll200rj ). 

Once a sink is formed, nearly all gas particles near the 
sinks will be accreted before approaching the d < r m i n cri- 
terion. However, secondary sinks can come within this dis- 
tance and become subject to strong N-body interactions. 
In physical reality, however, the sub-sink region will con- 
tain not only the protostar but also surrounding gas and an 
accretion disk. In addition, the protostars themselves can 
be very distended, with photospheres re aching ~ 1 AU (e.g. 
lOmukai fc Pallall2003l ; [Greif et alj|2012l : Smith et al. 2012a). 
Our accretion of all secondary sinks within d < r min there- 
fore allows us to account for sub-sink hydrodynamic affects 
which will mitigate those of pure N-body dynamics. 



2.5 Identifying Binaries 

Once a protostellar multiple system has formed, we iden- 
tify which pairs of sinks comprise a binary system. For each 
minihalo we determine which pairs of sinks are gravitation- 
ally bound by comparing their specific potential and kinetic 
energies, e p and e K . A pair of sinks is considered a binary 
system if e < 0, where e is the total orbital energy, 



e = e p + e k , 



-G(Mi+M 2 ) 



and 



1 2 

— V . 

2 



(2) 
(3) 



(4) 



pmax 



(i) 



Mi and M2 are the masses of the two sinks in question, r is 
the distance between the sinks, and v is the relative velocity 
between the sinks. Within each minihalo we find an average 
of between two and three binaries. 

Each binary is then assigned an overall mass Mbm = 
Mi + M2, a center-of-mass position, and a center-of-mass 
velocity. We next determine the semi-major axis, a, for each 
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Size: 7500 AU 




binary. This is derived by equating the above definition of e 
with the following energy equation: 

-GM hi n 



2a 

This yields 

_ -G(Mi+M 2 ) 



(->) 



(6) 



(7) 



2 (e k + e p ) 

We furthermore calculate the binary's mass ratio q, 
where 

_ Ma 
g ~~ Mi' 

In this case M2 corresponds specifically to the less massive 
of the binary members, and Mi corresponds to the more 
massive, such that the value of q is always between zero and 
one. 



3 RESULTS 

3.1 Initial Minihalo Collapse 

The initial collapse of gas to protostellar density exhibits 
substantial similarity within each minihalo. Figure [1] shows 
the thermal state of the star-forming gas within each mini- 
halo just prior to the formation of the first high-resolution 
sink, while Figure [2] shows their velocity profiles. At high 
densities of n > 10 s cm -3 , H2 three-body formation becomes 
rapid. This leads to nearly isothermal evolution as increased 
H2 cooling counteracts adiabatic heating while the gas con- 
denses to n > 10 12 cm -3 . The minihaloes show a striking 



similarity in thermal evolution in this density range. In the 
evolution of gas to the densities shown in Figure [T] (n > 
10 6 cm -3 ), we also find that the HD cooling is generally 
unimportant compared to H2 cooling, 

Figure [3] shows the angular momentum profile of the 
gas in the maximally resolved simulations just before the 
first sink has formed, measured with respect to the center- 
of-mass of all gas included in the 10 pc simulations. These 
profiles all have roughly the same magnitude and powerlaw 
slope, such that j(r) oc M enc , where M enc is the enclosed 
gas mass within a given radius r from the center-of-mass. 
This convergence is found not only among minihaloes in our 
simulation, but also including those from differing cosmolog- 
ical realizations, performed with varying m e thods of com- 
putat ional hydrodynamics (|Abel et al 1 120021 ; lYoshida et al.l 
l200rj ). These angular momentum distributions may addi- 
tionally lead to h i gh sp in rates of individual protostars (e.g. 
IStacv et al.ll201ll . l2012D . 

The density structure of the central gas can be seen 
in Figure [4] for a representative number of the minihaloes. 
The similarity between minihaloes is again apparent, though 
some of the star-forming gas exhibits greater elongation and 
a more disk-like structure. Similar variation is seen when we 
compare the velocity and turbulent structure in Fig. [21 We 
measure M tur b over a range of radial bins according to the 
following: 



,,2 2 mi 1 ^2 

M turb C s = -TT [Vi ~ Wrot - « rad ) , 



(8) 



where c s is the sound speed of the radial bin, mi is the mass 
of a gas particle contributing to the bin, and M is the total 
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Figure 5. History of the mass growth of the stellar systems in each minihalo. Solid lines represent the total sink mass. Dashed lines 
show the growth of the most massive sink. Dotted lines depict the growth of the second-most massive sink. Average accretion rates vary 
by over a factor of three between minihaloes. 



mass of the bin. Between 20 and 10,000 AU, both the rota- 
tional velocity v TO t and radial velocity i; ra d remain at a few 
km s _1 , while the turbulence remains approximately sonic, 
never reaching more than twice the thermal velocity of the 
gas. In each case, both radial and infall velocities are on the 
order of ~ 1/2 of the free-fall velocity vg. However, the rela- 
tive magnitude of i> ra d and v ro t varies, with w ra( j dominating 
in some cases and w ro t dominating in others. 

After the central gas in each minihalo reaches the den- 
sity threshold for sink formation, in all cases a stellar mul- 
tiple system quickly forms. This fra gmentation is in agree- 
ment with the parameter study of iMachida et all l|2008|) . 
They modeled their initial conditions as modified Bonnor- 
Ebert spheres, which well-describes collapsing and approxi- 
mat ely isothermal gas like that within star-forming regions 
(e.g.lEbertJll955l ; iBonnoJ 19561 : IStacv et al.ll2010l : iGreif et al.1 
and found the following fragmentation condition: 

Qo > n cri t = 4 x 10~ 17 ( - n , ?1 ° , V'V 1 , (9) 
VlO^cm^/ 

where f2o an d no are the angular velocity and central num- 
ber density of the initial core of the Bonnor-Ebert sphere. 
We determine $1q of the central core when no is between 



10 3 cm 3 and 10 4 cm 3 , corresponding to fi C rit values be- 
tween > 4 x 10~ 17 s" 1 and ~ 2 x 10" 16 s" 1 . We take a 
mass- weighted average of the angular velocity of all gas par- 
ticles with n > 10 3 cm -3 and find in each case that f2o 
ranges between two to three orders of magnitude greater 
than ficrit- The subsequent flattening and fragmentation of 
th e dense gas is thus con sistent with the criterion presented 
bv lMachida et al.1 (|200ST >. 



3.2 Mass Function 

The mass growth history of the sinks shows substantial vari- 
ation among the ten minihaloes, as is shown in Figures[5]and 
[13 At the end of the simulations, the total mass accreted by 
all sinks ranges from 30 — 120 Mq for a given minihalo (Fig. 
[5|- The total mass of gas within each minihalo is typically 
< 10 5 Mq, implying a star- formation efficienty of ~ 10 -3 . 
The distribution of sink masses after 5000 yr of accretion 
can be seen for the individual minihaloes in Figure |6j where 
the variation between minihaloes is even more striking. In 
particular, while some minihaloes have only one or two sinks 
with M* < 1 Mq, this number can be as large as nine. We 
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Figure 6. Distribution of sink masses in each minihalo after ~ 5000 yr of accretion. Dashed red lines depict powerlaw fits to each 
distribution, dN/dm oc M^ a . We find that a ranges from zero (Minihaloes 1 and 2) to 0.57 (Minihalo 10). Dashed blue lines show an 
example distribution for a = 2, the minimum powerlaw that will yield a top-heavy mass function. Each a = 2 distribution is normalized 
to have the same total mass as that found within the sinks of the corresponding minihalo. 



find a range of powerlaw fits to the relation AN /Am oc M^ a , 
where N is the number of sinks that lie within a logarith- 
mically scaled mass bin. This range extends from a = for 
Minihaloes 1, 2 and 9, to a = 0.57 for Minihalo 10. Consid- 
ering the distribution across all minihaloes (Fig. [7]), we find 
AN /Am oc Mr ' 17 . 

Unlike the observed present-day initial mass function 
(IMF) , which is characterized by a Salpeter slope of a = 2.35 
towards the high mass end, our simulations yield a top-heavy 
mass function in that the majority of the protostellar mass 
resides within the most massive protostars. Top-heaviness 
applies to mass functions with a < 2, since this is the re- 
quirement for the average stellar mass to be domina ted by 
the upper limit of the mass range (e.g. lBromrr]|2012l ). If we 
consider smaller values of a to be more top-heavy, then our 
fitted mass function is less top -heavy than that f ound after 
~ 1000 yr of sink accretion in iGreif et al.1 (|201ll ). who find 
a < 0. 

Including radiative feedback would likely modify our 
values for a, though the magnitude of the modification is 
unce rtain. While feedb ack does not prevent fragmentation 
(e.g. ISmith et al]|201ll . Stacy et al. 2012b), heating of the 



gas may suppress low-mass fragmentation to y ield smaller 
a valu es closer to zero (i.e., a flatter spectrum). IGreif et al.l 
l)201ll ) resolved smaller mass scales and again found smaller 
values of a. However, following such a resolved simulation 
for longer than their 1000 yr may reveal a population of very 
low-mass stars that arises at later times, instead increasing 
a. Our simulations do not resolve down to the opacity limit 
of fragmentation, the scale on which a small ~ 10~ 2 Mp> 
hydrostatic core first f orms (e.g. lReej[l976l : IOmukai fc Nishil 
ll99cj|G~reif et al1l2012l ). We may therefore underestimate the 
number of such low-mass protostellar fragments that form. 



3.3 Binary Fraction and Merger Rate 

Every minihalo contains at least one binary, and on average 
each halo hosts between two and three binaries. We define 
the binary fraction as 

B 



/b 



(10) 



S + B' 

where S and B are the number of single and binary sys- 
tems, respectively. Summing over all sinks remaining in each 
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Figure 7. Total distribution of sink masses in all minihalocs. 
Red dashed line depicts a powerlaw fit to the overall distribution, 
dN/dm oc M~ ' 17 . Blue dotted lines depict the total mass dis- 
tribution from the 'merging sink' simulations of Greif et al. 201f 
and a powcrlow fit to this distribution (a = —0.17). 



minihalo at the end of the simulations, we have S = 42 
and B = 24, yielding /b = 0.36. Calculated slightly dif- 
ferently, this means any individual star has a probability 
2B/(S + 2B)=53% of having a binary companion. This 
is a substantial fraction, and similar to the range of mas- 
sive star binary fractions observed in various clusters in the 
Milky Way fe.g. lSana fc Evans!l201ll ;l Kiminki fc Kobulnickvl 



Milky 
2012). 



While a total of 217 sinks form in our suite of mini- 
haloes, 90 remain after 5000 yr. Over than half (59%) of 
the sinks are therefore lost to mergers. T his is similar to 
the merger rate found in iGreif et all l|2012t ). However, it will 
require higher-resolution simulations to determine which of 
these mergers were actually unresolved but long-liv e d tigh t 
(a < 20 AU) binaries. The results of IGreif et all (|2012h . 
as well as our criterion that only non-rotationally supported 
sinks are merged, support the likelihood that a large number 
of these sink mergers were indeed true protostellar mergers. 



|2007| ; ISana fc Evans] l201ll ; | Kiminki fc Kobulnickvl 120121 ). 

However, the fit found from our simulation is significantly 
steeper. As will further be discussed in Section 4, there are 
multiple reasons for this. The spatial resolution of the simu- 
lations was limited, so the distributions cannot sample orbits 
with a < 20 AU. This contributes to the apparent lack of 
hard binaries. In addition, the binaries have only evolved 
for a few thousand years, and their orbits may evolve and 
harden over longer periods of time. 

The distribution of mass ratio q (Figure |9j goes as 
dN/dq oc g -0 ' 55 , while for the cumulative distribution we 
find a best-fit powerlaw of F q oc q ' 75 . The average value of 
q is 0.35. We comp are our q distribution to tha t found for 
solar-type stars by iDuquennov fc Mayor] (|l99ll ). shown as 
the blue dotted line in Fig. [9] They observed a q-distribution 
peaked at 0.23, while our distribution is peaked at q < 0.1, so 
our simulations have a comparatively high fraction of low-q 
binaries. 

Our rate of low-g binaries is also high when compared 
with m ore massive star cluster s in th e Milky Way. For in- 
stance, iKiminki fc Kobulnickvl |2012l ) find that dN/dq oc 
q ' 1 , while our fit falls off mu ch more steeply and i s slightly 
outside their range of error. I Sana fc Evans! (|201ll ) present 
observations which yield a uniform distribution of q down 
to q — 0.2. This corresponds to the cumulative distribution 
shown as the blue line in the right panel of Fig. [9] 

If the sink particles of the simulations continued to ac- 
crete mass for longer periods of time, it is possible that 
the q distribution would shift toward a larger percentage 
of binaries with nearl y equal mass. As described in, e.g., 
iBate fc Bonnelj (| 19971 ). if gas with sufficient angular mo- 
mentum flows onto a binary star, it will accrete preferen- 
tially onto the less massive companion. Angular momentum 
conservation prevents rotating gas from reaching the larger 
star without the aid of, e.g., gravitational torques. The gas 
instead remains at some radius beyond the more massive 
star where it is more easily captured by the less massive 
companion, resulting in a tendency for q to evolve towards 
a value of one. This is similar to the fragm entation-induced 
starvation described bv lPeters et alj (|2010l ). In their simula- 
tions of star-forming molecular clouds, they found that gas 
inflowing through a stellar disk was accreted by lower-mass 
stars before it could reach the more massive and centrally- 
located star. 



3.4 Orbital Period and Mass Ratio Distribution 

The characteristics of the binary distribution are shown 
in Figures [8] and [9] The distribution of the binary pe- 
riods is peaked at 10 5,J days, or 870 yr. In compari- 
son with nearby solar-type stars (blue dotted line in Fig. 
03, this peak is shifted to longer periods. For instance, 
(fpucmennov fc Mavor|[T99ll ) observed a P distribution that 
is peaked at 10 4 ' 8 days, or 170 years. This corresponds to 
a distribution of the semi-major axis, a, that is peaked 
around ~ 250 AU, with an average value of 760 AU. In 
Figure [8] we also show the related cumulative distribution 
of the binary periods, Fp, normalized such that Fp(P < 
00) = 1. We display a best-fit powerlaw for this distrib ution, 
fp oc log(P) 2 ' 54 , and compare to a 'Opik distribution' (Opik; 
Il925l ). where Fp oc log(P). Opik's law i s consistent with ob- 
servations of massive star binaries (e.g. iKouwenhoven et al.l 



3.5 Protostellar Ejection Rate 

Of the total of 90 sinks remaining at the end of our simu- 
lation, 39 have radial velocities that are greater than their 
escape velocities (Figure [TO)) , where we define v BS c as 



GM c: 



(11) 



where M Gnc is the mass enclosed between the center of mass 
(COM) and the sink particle at distance r from the COM. 
The COM is determined using all gas and sink particles 
with n > 10 10 cm -3 . Roughly half (43%) of the sinks can 



thus leave t he stellar disk 
IGreif et al.1 



similar to the rate found in, e.g. 



20 111) in th e ir sim ulations of primordial star for- 



Bate et al.1 (|2003h in their numerical studies of 



mation and 
present-day, low-mass star formation. 
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Figure 8. Left: Distribution of the log of the period, P, for the binaries found within each minihalo. Note that units of P is in days. Red 
dashed line shows a Gaussian fit to the distribution. Blue dotted line shows an example fit to data for solar-type stars (Duqucnnoy & 
Mayor 1991). The distribution of the binary periods is peaked at 10 5,5 days, or 870 yr. Right: Cumulative distribution Fp of the orbital 
period of the binaries identified in all minihaloes. Red dashed line shows a powerlaw fit, Fp oc log(P) 2,54 , while the blue line shows a 
sample Opik distribution for comparison. 




Figure 9. Left: Distribution of the mass ratio q between the members of binaries located in each minihalo. Red line is the powerlaw fit 
dN/dq oc q~ °- 55 . Blue dotted line shows an example fit to data for solar-type stars (Duqucnnoy & Mayor 1991). The average value of q is 
0.34. Right: Cumulative distribution F q of the mass ratio of all binaries in our simulation. Red dashed line is a powerlaw fit, F q oc q 075 . 
The blue line depicts a uniform distribution down to q = 0.2, as observed in massive star clusters in the Milky Way (e.g. Sana & Evans 
2011). 



In addition, we define a second escape velocity, that 
required to exit the minihalo: 



Wesc.haio = \ / ~5kms , (12) 

V ''halo 

where we define Mhaio = 4 x 10 5 Mq and rhaio = 80 pc, 
the average values for the minihaloes of our suite. Whether 
v csc or Ueschaio is larger varies depending on the sink and its 
location. A total of 53 sinks have u r ad > iw.haio, which is 
greater than the 39 which can escape from the stellar disk. 
Thus, the gravity of the stellar disk provides the stronger 
escape condition. In the majority of cases, for a sink to have 
sufficient velocity to escape the disk, it must already have 
sufficient velocity to escape the minihalo (left panel of Fig- 



ure [TO}. For a sink traveling at tw.haio, the time to cross 
the virial radius of the minihalo is ~ 2 x 10 7 yr. Note that 
within this time the minihalo is not expected to undergo sig- 
nificant mass growth it self, no more than d oubling in mass 
(e.g. Stacy et al. 2011b. IWechsler et al.ll2002T ). Thus, iw.haio 
will not increase by more than 1-2 km s _1 . 

Of the ejected protostars, 18 are of mass M„ < 1 Mq 
and thus have a possibility of surviving to the present day 
(Clark et al 2011b). This depends, however, on whether the 
stars will accrete more mass at a later time and thus become 
shorter lived. Identifying them as Pop III stars would require 
that, as they become incorporated into the Milky Way or 
nearby dwarf galaxies, they will not accrete any enriched ma- 
terial that will remain in their atmospheres and mask them 
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Figure 10. Left: Ratio of |u ra d| to t> eS c for each sink particle. Sinks with M* < 1 Mq are denoted with filled circles, while diamonds 
represent sinks with M* > 1 Mq. Horizontal dashed line shows where \v I!id /v eEC \ = 1. Sinks above the horizontal dashed line are able to 
escape from the stellar disk. The vertical dashed line denotes v eac halo> the velocity necessary to escape the minihalo. Right: Variation of 
l^radl with radius. Dashed line again shows « C sc,halo- 



as Po p II stars (e.g. iFrebel et al.ll2009l : Ijohnson fc KhochfaJ enhance the physical realism of subseqent simulations by 
l201lf ). including these processes and resolving smaller scales. 



4 CAVEATS 

The results presented in the simulation have a number of 
caveats that should be noted. First, these simulations fol- 
lowed only 5000 yr of protostellar mass accretion. The dis- 
tribution of sink masses and binary properties will continue 
to evolve beyond this point, as some of the sinks will con- 
tinue to grow in mass and orbits may harden or undergo 
disruption through close encounters with other sinks. We 
also cannot study binaries within the resolution limit of the 
simulation, 20 AU. Thus, the sink mergers discussed previ- 
ously may have instead become very tight and unresolved 
binaries. Objects formed from subsequent fragmentation on 
sub-sink scales also could not be addressed in this study. If 
we did presume that all sink mergers instead became tight 
binaries orbiting with semi-major axis of 10 AU, we would 
have a binary fraction of 73%. The average orbital radius 
of all binaries would then be 170 AU, and the average mass 
ratio would still be 0.34. 

In addition, we do not include protostellar feedback 
in our simulations. The Hb-dissociating and ionizing radi- 
ation emitted from the growing protostars may substan- 
tially change the characteristics of the stellar disk and 
eventually cut off mass flow onto the d isk altogether (e.g. 
iHosokawa et alJbOlllJSmith et alj|201ll ; Stacy et al. 2012b). 
Another important physical process not included in our sim- 
ulations was magnetic fields, which may affect the angular 
momentum build-up of the protostellar disk, the subsequent 
fragmentation, and the mass flow onto the protostars (e.g 
Tan fc Blackmanl 12004 ISilk fc LangeJ 12009; iMaki fc Susa 



20071 : ISchleicher et al.ll2010l ; lFederrath et al.ll201ll ; ISur et a! 
20121 ; Schober et al. 2012a, 2012b). In the future we plan to 



5 SUMMARY AND CONCLUSIONS 

We have examined the formation of Pop III multiple sys- 
tems within ten different minihaloes extracted from a 1.4 
Mpc (comoving) cosmological box. The properties of each 
stellar group varies between minihaloes, with mass functions 
ranging from flat to a = 0.57, and average stellar accretion 
rates varying by over a factor of three between minihaloes. 

Looking at the combined properties of the ten stellar 
systems, we find that 50% of the protostars merge with 
larger protostars, while 50% of the remaining protostars are 
ejected from their stellar disk. We furthermore find a bi- 
nary fraction of ~ 36%, translating to a 53% probability 
that any given star will have a companion. Some of these 
binaries have semi-major axes that extend to 3000 AU. The 
distribution of semi-major axes peaks slightly at ~ 300 AU, 
while the mass ratio is nearly uniformly distributed between 
zero and one. Recent work by Smith et al. (2012b) implies 
that this result may be modified by the effects of heating 
and ionization from DM annihilation (DMA). If their find- 
ings hold up to further investigation, DMA may suppress 
fragmentation and lead to wider binaries of a > 1000 AU. 

Our simulation did not have sufficient resolution to re- 
solve protostellar mergers or the formation of tight binaries, 
as these processes occur on scales smaller than our resolution 
length of 20 AU. If a significant fraction of these sink merg- 
ers did represent true protostellar mergers or the formation 
of contact binaries, this has imp ortant implic ations for the 
Pop III stellar evolution (see, e.g. lLangei)2012l and references 
therein). During a common-envelope phase, H-envelope re- 
moval by a close binary companion may allow for a star as 
low as 10 Mp) to beco me a WR star (e.g. IVanbeveren et all 
1 19981 ; ICrowthedl2007h . A stellar merger may rejuvenate the 
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mass-gai ning star, making it appear younger than its true 
age ( e.g. lHellingslll983l ; H rami fc Langerlll995l ; lDrav fe Tout! 
120071 ). Alternativ ely, a star with an un dermassive helium 
core may result (|Braun fc Langerl Il995l ). Such stars may 
avoid the red supergiant phase and instead explode as 
blue supergiants. Mass accretion through Roche-lobe over- 
flow may also spin up a star, allowing for rotationally- 
induced mixing and possibly chemically homogeneous evo- 
lution, thereby altering th e star's luminosity, temperature, 
and m etal production (e.g. lEkstrom et al.ll2008l ; lYoon et al.l 
2012). In particular, rotationally induced mixing may al- 
ter the nucleosynthesis within a within a low-metallicity or 
primordial massive star to yield greater CNO-element abun- 
danc es, especially 14 N (e.g. lEkstrom et al1l2008l ; lYoon et al.l 
120121 '). This additionally leads to a neutron excess in the 
core, enhancing the production of s-process elements (e.g. 
IPignatari et al.ll2008ft . 

The binarity of Pop III stars also provides an important 
pathway for the formation of GRBs and high-mass X-ray bi- 
naries. If one of the binary members is in the mass range to 
collapse into a BH, then a close companion may remove its 
hydrogen envelope and allow for penetration of a G RB jet 
to the surface (e.g. lLee et~ai]|2002l ; llzzard et al.ll2004l '). Once 
the BH remnant is left behind, the Roche lobe overflow from 
the companion enables the BH to sustain a high luminosity 
and X-ray flux. These sources could significantly contribute 
to reionization, and exert feedback effects during the as- 
semb ly of the first galaxies (|Mirabel et alj|201ll ; Ijeon et "all 
I2012D . 

Pop III binarity is thus a crucial aspect to understand- 
ing Pop III evolution and feedback. A significant percent- 
age of Pop III stars exist in binaries and undergo stellar 
mergers. Greater computational power will allow for future 
simulations which include more physical processes such as 
feedback and magnetic field growth, and which will also fol- 
low the evolution of Pop III multiple systems with increased 
resolution over longer periods of time. Combined with ob- 
servational constraints, this will yield an increasingly refined 
understanding of Pop III stars and their role in the early 
Universe. 



ACKNOWLEDGMENTS 

The authors wish to thank John Mather for insightful dis- 
cussion during the development of this work. AS is grateful 
for support from the JWST Postdoctoral Fellowship through 
the NASA Postdoctoral Program (NPP). VB acknowledges 
support from NASA through Astrophysics Theory and Fun- 
damental Physics Program grant NNX09AJ33G and from 
NSF through grant AST-1009928. Resources supporting this 
work were provided by the NASA High-End Computing 
(HEC) Program through the NASA Advanced Supercom- 
puting (NAS) Division at Ames Research Center. 



REFERENCES 

Abel T., Bryan G. L., Norman M. L., 2002, Sci, 295, 93 
Alvarez M. A., Bromm V., Shapiro P. R., 2006, ApJ, 639, 
621 

Alvarez M. A., Wise J. H., Abel T., 2009, ApJ, 701, L133 



Bate M. R., Bonnell I. A., 1997, MNRAS, 285, 33 
Bate M. R., Bonnell I. A., Bromm V., 2003, MNRAS, 339, 
577 

Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 
362 

Bonnor W. B., 1956, MNRAS, 116, 351 
Braun H., Langer N., 1995, A&A, 297, 483 
Bromm V., 2012, Second Workshop on Numerical and Ob- 
servational Astrophysics: From the First Structures to the 
Universe Today (arXiv: 1203.3824) 
Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23 
Bromm V., Kudritzki R. P., Loeb A., 2001, ApJ, 552, 464 
Bromm V., Loeb A., 2003, ApJ, 596, 34 
Bromm V., Loeb A., 2004, New Astron., 9, 353 
Bromm V., Yoshida N., Hernquist L., 2003, ApJ, 596, L135 
Brott I., de Mink S. E., Cantiello M., Langer N., de Koter 
A., Evans C. J., Hunter I., Trundle C, Vink J. S., 2011, 
A&A, 530, A115 
Cantiello M., Yoon S.-C, Langer N., Livio M., 2007, A&A, 
465, L29 

Chatzopoulos E., Wheeler J. C, 2012, ApJ, 748, 42 
Clark P. C, Glover S. C. O., Klessen R. S., 2008, ApJ, 672, 
757 

Clark P. C, Glover S. C. O., Klessen R. S., Bromm V., 

2011a, ApJ, 727, 110 
Clark P. C, Glover S. C. O., Smith R. J., Greif T. H., 

Klessen R. S., Bromm V., 2011b, Sci, 331, 1040 
Crowther P. A., 2007, ARA&A, 45, 177 
de Mink S. E., Cantiello M., Langer N., Pols O. R., Brott 

I., Yoon S.-C, 2009, A&A, 497, 243 
Donati J.-F., Landstreet J. D., 2009, ARA&A, 47, 333 
Dray L. M., Tout C. A., 2007, MNRAS, 376, 61 
Duquennoy A., Mayor M., 1991, A&A, 248, 485 
Ebert R., 1955, Zeitschrift ft Astrophysik, 36, 222 
Ekstrom S., Meynet G., Chiappini C, Hirschi R., Maeder 

A., 2008, A&A, 489, 685 
Federrath C, Chabrier G., Schober J., Banerjee R., Klessen 

R. S., Schleicher D. R. G., 2011, Physical Review Letters, 

107, 114504 

Ferrario L., Pringle J. E., Tout C. A., Wickramasinghe 

D. T., 2009, MNRAS, 400, L71 
Frebel A., Johnson J. L., Bromm V., 2009, MNRAS, 392, 

L50 

Greif T., Springel V., White S., Glover S., Clark P., Smith 
R., Klessen R., Bromm V., 2011, ApJ, 737, 75 

Greif T. H., Bromm V., Clark P. C, Glover S. C. O., Smith 
R. J., Klessen R. S., Yoshida N., Springel V., 2012, MN- 
RAS, 424, 399 

Greif T. H., Glover S. C. O., Bromm V., Klessen R. S., 

2010, ApJ, 716, 510 
Greif T. H., Johnson J. L., Bromm V., Klessen R. S., 2007, 

ApJ, 670, 1 

Greif T. H., Johnson J. L., Klessen R. S., Bromm V., 2009, 

MNRAS, 399, 639 
Haiman Z., 2011, Nat, 472, 47 

Haiman Z., Thoul A. A., Loeb A., 1996, ApJ, 464, 523 
Hamann W.-R., Grafener G., Liermann A., 2006, A&A, 
457, 1015 

Heger A., Woosley S. E., 2002, ApJ, 567, 532 
Hellings P., 1983, Ap&SS, 96, 37 

Hosokawa T., Omukai K., Yoshida N., Yorke H. W., 2011, 
Science, 334, 1250 



© 2012 RAS, MNRAS 000.ITHT41 



Pop III Binaries 13 



Izzard R. G., Ramirez-Ruiz E., Tout C. A., 2004, MNRAS, 
348, 1215 

Jeon M., Pawlik A. H., Greif T. H., Glover S. C. O., Bromm 

V., Milosavljevic M., Klessen R. S., 2012, ApJ, 754, 34 
Johnson J. L., Bromm V., 2007, MNRAS, 374, 1557 
Johnson J. L., Greif T. H., Bromm V., 2007, ApJ, 665, 85 
Johnson J. L., Khochfar S., 2011, MNRAS, 413, 1184 
Karlsson T., Bromm V., Bland-Hawthorn J., 2012, Rev. 

Modern Phys., in press (arXiv: 1101.4024) 
Kiminki D. C, Kobulnicky H. A., 2012, ApJ, 751, 4 
Kitayama T., Yoshida N., Susa H., Umemura M., 2004, 
ApJ, 613, 631 

Kouwenhoven M. B. N., Brown A. G. A., Portegies Zwart 

S. F., Kaper L., 2007, A&A, 474, 77 
Kowalska I., Bulik T., Belczynski K., 2012, A&A, 541, A120 
Kuhlen M., Madau P., 2005, MNRAS, 363, 1069 
Langer N., 2012, ARA&A, 50, 107 

Langer N., Cantiello M., Yoon S.-C, Hunter I., Brott I., 
Lennon D., de Mink S., Verheijdt M., 2008, in Bresolin F., 
Crowther P. A., Puis J., eds, IAU Symposium Vol. 250, 
Rotation and Massive Close Binary Evolution, p. 167 

Lee C, Brown G. E., Wijers R. A. M. J., 2002, ApJ, 575, 
996 

Lee W. H., Ramirez-Ruiz E., 2006, ApJ, 641, 961 
Linden T., Kalogera V., Sepinsky J. F., Prestwich A., Zezas 

A., Gallagher J. S., 2010, ApJ, 725, 1984 
Machida M. N., Omukai K., Matsumoto T., Inutsuka S.-i., 

2008, ApJ, 677, 813 
Madau P., Ferrara A., Rees M. J., 2001, ApJ, 555, 92 
Madau P., Rees M. J., Volonteri M., Haardt F., Oh S. P., 

2004, ApJ, 604, 484 
Maeder A., Meynet G., 2012, Rev. Modern Phys., 84, 25 
Maio U., Khochfar S., Johnson J. L., Ciardi B., 2011, MN- 
RAS, 414, 1145 
Maki H., Susa H., 2007, PASJ, 59, 787 
Martel H., Evans N. J., Shapiro P. R., 2006, ApJS, 163, 
122 

Mason B. D., Hartkopf W. I., Gies D. R., Henry T. J., 

Helsel J. W., 2009, AJ, 137, 3358 
Meynet G., Maeder A., 2005, A&A, 429, 581 
Milosavljevic M., Bromm V., Couch S. M., Oh S. P., 2009, 

ApJ, 698, 766 

Mirabel I. F., Dijkstra M., Laurent P., Loeb A., Pritchard 

J. R., 2011, A&A, 528, A149 
Mori M., Ferrara A., Madau P., 2002, ApJ, 571, 40 
Nakauchi D., Suwa Y., Sakamoto T., Kashiyama K., Naka- 

mura T., 2012, ApJ, in press (arXiv: 1207.2835) 
Norman M. L., O'Shea B. W., Paschos P., 2004, ApJ, 601, 

L115 

Omukai K., Nishi R., 1998, ApJ, 508, 141 

Omukai K., Palla F., 2003, ApJ, 589, 677 

Opik E., 1925, Tartu Obs. Publ., 26, 1 

Palla F., Salpeter E. E., Stahler S. W., 1983, A[J, 271, 632 

Peters T., Banerjee R., Klessen R. S., Mac Low M.-M., 

Galvan-Madrid R., Keto E. R., 2010, ApJ, 711, 1017 
Petrovic J., Langer N., van der Hucht K. A., 2005, A&A, 

435, 1013 

Petrovic J., Langer N., Yoon S.-C, Heger A., 2005b, A&A, 
435, 247 

Pignatari M., Gallino R., Meynet G., Hirschi R., Herwig 

F., Wiescher M., 2008, ApJ, 687, L95 
Rees M. J., 1976, MNRAS, 176, 483 



Ricotti M., Ostriker J. P., 2004, MNRAS, 352, 547 

Sana H., de Mink S. E., de Koter A., Langer N., Evans 
C. J., Gieles M., Gosset E., Izzard R. G., Le Bouquin 
J.-B., Schneider F. R. N., 2012, Sci, 337, 444 

Sana H., Evans C. J., 2011, in Neiner C, Wade G., Meynet 
G., Peters G., eds, IAU Symposium Vol. 272 of IAU Sym- 
posium, The multiplicity of massive stars, pp 474-485 

Schaerer D., 2002, A&A, 382, 28 

Schleicher D. R. G., Banerjee R., Sur S., Arshakian T. G., 
Klessen R. S., Beck R., Spaans M., 2010, A&A, accepted 
(arXiv: 1003.1135) 

Schober J., Schleicher D., Federrath C, Glover S., Klessen 
R. S., Banerjee R., 2012a, ApJ, 754, 99 

Schober J., Schleicher D., Federrath C, Klessen R., Baner- 
jee R., 2012b, PRE, 85, 026303 

Silk J., Langer M., 2006, MNRAS, 371, 444 

Smith R. J., Glover S. C. O., Clark P. C, Greif T., Klessen 
R. S., 2011, MNRAS, 414, 3633 

Smith R. J., Hosokawa T., Omukai K., Glover S. C. O., 
Klessen R. S., 2012a, MNRAS, 424, 457 

Smith R. J., Iocco F., Glover S. C. O., Schleicher D. R. C, 
Klessen R. S., Hirano S., Yoshida N., 2012b, ApJ, submit- 
ted (arXiv: 1210.1582) 

Sokasian A., Yoshida N., Abel T., Hernquist L., Springel 
V., 2004, MNRAS, 350, 47 

Springel V., 2005, MNRAS, 364, 1105 

Stacy A., Bromm V., Loeb A., 2011, MNRAS, 413, 543 

Stacy A., Bromm V., Loeb A., 2011b, ApJ, 730, LI 

Stacy A., Greif T. H., Bromm V., 2010, MNRAS, 403, 45 

Stacy A., Greif T. H., Bromm V., 2012b, MNRAS, 422, 290 

Stacy A., Greif T. H., Klessen R. S., Bromm V., Loeb A., 
2012, MNRAS, submitted (arXiv: 1209. 1439) 

Sur S., Federrath C, Schleicher D. R. G., Banerjee R., 
Klessen R. S., 2012, MNRAS, 423, 3148 

Suwa Y., Ioka K., 2011, ApJ, 726, 107 

Swartz D. A., Soria R., Tennant A. F., 2008, ApJ, 684, 282 

Tan J. C, Blackman E. G., 2004, ApJ, 603, 401 

Tegmark M., Silk J., Rees M. J., Blanchard A., Abel T., 
Palla F., 1997, ApJ, 474, 1 

Tohline J. E., 2002, ARA&A, 40, 349 

Tornatore L., Ferrara A., Schneider R., 2007, MNRAS, 382, 
945 

Turk M. J., Abel T., O'Shea B., 2009, Sci, 325, 601 

Turk M. J., Clark P., Glover S. C. O., Greif T. H., Abel 

T., Klessen R., Bromm V., 2011, ApJ, 726, 55 
Vanbeveren D., De Loore C, Van Rensbergen W., 1998, 

A&ARv, 9, 63 
Wada K., Venkatesan A., 2003, ApJ, 591, 38 
Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov 

A. V., Dekel A., 2002, ApJ, 568, 52 
Wellstein S., Langer N., Braun H., 2001, A&A, 369, 939 
Whalen D., Abel T., Norman M. L., 2004, ApJ, 610, 14 
Wise J. H., Abel T., 2008, ApJ, 685, 40 
Woosley S. E., 1993, ApJ, 405, 273 
Woosley S. E., Heger A., 2006, ApJ, 637, 914 
Yoon S., Langer N., 2005, A&A, 443, 643 
Yoon S.-C, Dierks A., Langer N., 2012, in press (arXiv: 

201.2364) 

Yoshida N., Abel T., Hernquist L., Sugiyama N., 2003, 
ApJ, 592, 645 

Yoshida N., Omukai K., Hernquist L., 2008, Sci, 321, 669 



© 2012 RAS, MNRAS OOO.QjITj] 



14 A. Stacy and V. Bromm 

Yoshida N., Omukai K., Hernquist L., Abel T., 2006, ApJ, 
652, 6 

Zhang W., Woosley S. E., Heger A., 2004, ApJ, 608, 365 



© 2012 RAS, MNRAS OOO.ITHT41 



